---
title: "Oficina DiD"
author: "Patrícia Pain"
date: "`r Sys.Date()`"
output: html_document
---
# Instalar e carregar pacotes

```{r}
pacotes <- c("basictabler","car","caret","conflicted","correlation","cowplot","dplyr","faraway","fastDummies","flextable", "foreign","ggrepel","ggtree","ggplot2","ggpubr","graphics", "grid", "gridExtra","gtsummary","Hmisc","jsmodule", "jtools","knitr","kableExtra","lmtest", "lubridate","magick", "MASS","mfx","minqa","mgcv","nnet","plotly", "palmerpenguins","performance","pglm",  "pgls","plm","pROC","pscl","psych","readr","readxl","ReporteRs", "rgl","reshape2", "ROCR","scales","sjlabelled","stargazer","stats","tidyr","tidyverse","tinytex",  "tseries","visreg","xlsx","xtable","writexl","vdr","MatchIt")
if(sum(as.numeric(!pacotes %in% installed.packages())) != 0){
  instalador <- pacotes[!pacotes %in% installed.packages()]
  for(i in 1:length(instalador)) {
    install.packages(instalador, dependencies = T)
    break()}
  sapply(pacotes, require, character = T)
} else {
  sapply(pacotes, require, character = T)
}
```

# Carregando Base

```{r}
load("C:/Users/Usuário/Desktop/Oficina DiD/Environmen inicial.RData")
```

Criando variáveis de evento e grupo de tratamento

# Variáveis utilizadas

### TRATAMENTO
País

Grupo de tratamento: empresas listadas nos EUA
Grupo de controle: empresas listadas no Japão

* Indicadores de ambiente institucional com características semelhantes
Developed Countries
Common Law
High Accounting Enforcement

### EVENTO
CAM - Critical Audit Matters

### CONTROLES:
Big4
AUD
SIZE
ROA
LEVER
BTM
CFO
LCS
COV

```{r}
Base_DiD <- Base_DiD %>%
  dplyr::mutate(Evento = ifelse(Year>2018, 1,0),
                Pais = ifelse(Country == "United States of America", 1,0))

# Varificando NAs
apply(apply(Base_DiD,2,is.na),2,sum)

# Olhando para a base
summary(Base_DiD)
```

# Estatística descritiva

Comparando as variáveis entre os grupos de tratamento e controle.

```{r}
Base_DiD %>%
  dplyr::select(DA_DSS_w1.bc, SIZE_w1, ROA_w1, LEVER_w1, BTM_w1, CFO_control_w1, SP, Pais, CommonLaw, HighEnforcement, Evento, Big4, AUD, LCS, Emerg, COV) %>%
  dplyr::mutate(Pais = factor(Pais, labels = c(0, 1))) %>%
  gtsummary::tbl_summary(by = Pais,
                         statistic = list(all_continuous() ~ "{median}, {mean} ({sd})",
                                          all_categorical() ~ "{n} ({p}%)"),
                         digits = all_continuous() ~ 2,
                         missing = "no", 
                         missing_text = "Missing",
                                      ) %>%
  add_n() %>% 
  add_p(list(all_continuous() ~ "t.test",
             all_categorical() ~ "kruskal.test")) %>%
  modify_header(label ~ "**Variáveis**") %>%
  modify_caption("Tabela X. Estatística Descritiva") %>%
  italicize_levels() %>%
  bold_p() %>%
  as_flex_table() 
```

# Accruals Discricionários - Dechow (DA_DSS)

Criar bases das variáveis dependentes

```{r}
bdDA <- Base_DiD %>%
  dplyr::filter(DA_DSS_w1.bc != "NA")
```

## PSM para DA_DSS

### Estimando as empresas comparáveis

Identificando as empresas da amostra com maior probabilidade de serem comparáveis

```{r}
PSM_DAPais <- glm(Pais ~ SIZE_w1 + ROA_w1 + LEVER_w1 + BTM_w1 + CFO_control_w1, family = binomial(), 
               data = bdDA)

summary(PSM_DAPais)
```

Estimando as empresas comparáveis

```{r}
Predict_PSM_DA <- data.frame(pr_score = predict(PSM_DAPais, type = "response"),
                     "EUA" = PSM_DAPais$model$Pais)
head(Predict_PSM_DA)
```

Visualizando os grupos previstos

```{r}
#Visualizando os grupos previstos em dois quadros
Grafico_Distr_original <- paste("Distribuição dos grupos:", c("EUA", "JAPAO"))
Predict_PSM_DA %>%
  dplyr::mutate(EUA = ifelse(EUA == 1, 1, 0)) %>%
  ggplot(aes(x = pr_score)) +
  geom_histogram(color = "white") +
  facet_wrap(~EUA) +
  xlab("Score PSM") +
  theme_bw()

#Visualizando os grupos previstos no mesmo quadro
Grafico_Distr_original <- paste("Distribuição dos grupos:", c("EUA", "JAPAO"))
Predict_PSM_DA %>%
  dplyr::mutate(EUA = ifelse(EUA == 1, "EUA", "JAPAO")) %>%
  ggplot(aes(x = pr_score, fill = EUA)) +  # Use 'fill' para cores diferentes
  geom_histogram(color = "white") +
  xlab("Score PSM") +
  theme_bw() +
  labs(title = Grafico_Distr_original) +
  scale_fill_manual(values = c("EUA" = "blue", "JAPAO" = "red"))  # Defina cores personalizadas
```

### Usando a função "matchit" para criar a nova base

```{r}
mod_PSM_DAPais <- matchit(Pais ~ SIZE_w1 + ROA_w1 + LEVER_w1 + BTM_w1 + CFO_control_w1,
  method = "nearest", data = bdDA)

bdDAPais <- match.data(mod_PSM_DAPais)

dim(bdDAPais)
```

## Cálculo das difereças

```{r}
bdDA_EUA_PreAdocao <- bdDAPais %>%
  dplyr::filter(Pais==1 & Evento==0)

bdDA_EUA_PosAdocao <- bdDAPais %>%
  dplyr::filter(Pais==1 & Evento==1)

bdDA_JAPAO_PreAdocao <- bdDAPais %>%
  dplyr::filter(Pais==0 & Evento==0)

bdDA_JAPAO_PosAdocao <- bdDAPais %>%
  dplyr::filter(Pais==0 & Evento==1)
```

Média dos grupos

```{r}
mean(bdDA_EUA_PreAdocao$DA_DSS_w1.bc) #Tratamento-Pre Evento
mean(bdDA_EUA_PosAdocao$DA_DSS_w1.bc) #Tratamento-Pos Evento
mean(bdDA_JAPAO_PreAdocao$DA_DSS_w1.bc) #Controle-Pre Evento
mean(bdDA_JAPAO_PosAdocao$DA_DSS_w1.bc)#Controle-Pos Evento
```

Calculando a diferença entre os grupos

```{r}
DifTrat <- (mean(bdDA_EUA_PosAdocao$DA_DSS_w1.bc) - mean(bdDA_EUA_PreAdocao$DA_DSS_w1.bc))

DifContr <-(mean(bdDA_JAPAO_PosAdocao$DA_DSS_w1.bc) - mean(bdDA_JAPAO_PreAdocao$DA_DSS_w1.bc))

DifTrat - DifContr
```

Outra forma

```{r}
((mean(bdDA_EUA_PosAdocao$DA_DSS_w1.bc) - mean(bdDA_EUA_PreAdocao$DA_DSS_w1.bc)) - 
    (mean(bdDA_JAPAO_PosAdocao$DA_DSS_w1.bc) - mean(bdDA_JAPAO_PreAdocao$DA_DSS_w1.bc)))
```

## Regressão para DA_DSS

### Rodando o pooled 

```{r}
mod_DSSPais_pooled <- plm(formula = DA_DSS_w1.bc ~ factor(Pais)*factor(Evento) + SIZE_w1 + ROA_w1 + LEVER_w1 + BTM_w1 + CFO_control_w1 + factor(AUD) + factor(LCS) + factor(COV),
                      data = bdDAPais, 
                      index = c ("Ticker", "Year"), 
                      model = "pooling")

summary(mod_DSSPais_pooled)
```

### Rodando o Efeitos-Aleatórios 

```{r}
mod_DSSPais_EA <- plm(formula = DA_DSS_w1.bc ~ factor(Pais)*factor(Evento) + SIZE_w1 + ROA_w1 + LEVER_w1 + BTM_w1 + CFO_control_w1 + factor(AUD) + factor(LCS) + factor(COV),
                      data = bdDAPais, 
                      index = c ("Ticker", "Year"), 
                  model = "random")

summary(mod_DSSPais_EA)
```

### Rodando o Efeitos-Fixos 

```{r}
mod_DSSPais_EF <- plm(formula = DA_DSS_w1.bc ~ factor(Pais)*factor(Evento) + SIZE_w1 + ROA_w1 + LEVER_w1 + BTM_w1 + CFO_control_w1 + factor(AUD) + factor(LCS) + factor(COV),
                      data = bdDAPais, 
                      index = c ("Ticker", "Year"), 
                  model = "between")

summary(mod_DSSPais_EF)
```

### Escolhendo o Melhor Modelo

Pooled x EF

Teste F de Chow \*A hipótese nula é de que há igualdade nos interceptos
e nas inclinações para todos os indivíduos, caracterizando o modelo de
dados agrupados (pooled). Se valor p\<0,05, o modelo de Efeitos Fixos é
melhor do que o modelo Pooled.

```{r}
pFtest(mod_DSSPais_EF,mod_DSSPais_pooled) #EF
```

Pooled x EA

Teste de Breusch e Pagan \*A aceitação da hipótese nula implica que o
modelo de dados agrupados (pooled) é preferível. Se valor p\<0,05, o
modelo de Efeitos Aleatórios é superior ao modelo Pooled.

```{r}
plmtest(mod_DSSPais_pooled, type="bp") #EA
```

EF x EA

Teste de Hausmann \*Se o teste rejeitar a hipótese nula, o modelo de
Efeitos Fixos é o mais adequado. Se valor p\<0,05 o modelo de Efeitos
Fixos foi considerado superior ao modelo de Efeitos Aleatórios.

```{r}
phtest(mod_DSSPais_EF,mod_DSSPais_EA) #EA
```

### Analisando os pressupostos para o modelo escolhido 

Normalidade dos resíduos 

(Anderson-Darling) normality test H0: Há distribuição normal dos
resíduos - se p-value \< 0,05, rejeita-se a hipótese nula, logo não há
distribuição normal dos resíduos.

```{r}
ad.test(mod_DSSPais_EA$residuals)
```

Homocedasticidade dos resíduos 

(Breusch-Pagan) \*hipótese nula
(p-valor\>0,05) é a de que não há homocedasticidade nos resíduos

```{r}
bptest(mod_DSSPais_EA)
```

Multicolinearidade

Um VIF entre 5 e 10 indica alta correlação, o que pode ser problemático.
E se o VIF for acima de 10, você pode assumir que os coeficientes de
regressão estão mal estimados devido à multicolinearidade.

```{r}
car::vif(mod_DSSPais_EA)
```

Correlação serial 

(teste Breusch-Godfrey/Wooldridge) \*hipótese nula
(p-valor\>0,05), não há problemas de correlação serial nos dados.

```{r}
pbgtest(mod_DSSPais_EA) 
```

### Exportando os resultados com erros padrões robustos clusterizados

```{r}
export_summs(mod_DSSPais_EA, scala = F,
             model.names = c("EQ.1-DSS"),
             error_pos = c("same", "below", "right"),
             bold_signif = 0.1,
             borders = 2,
             outer_borders = 2,
             statistics = c(N = "nobs.1", R2 = "r.squared",
                            adj.R2 = "adj.r.squared", p.value = "p.value", 
                            "GL" = "df", AIC = "AIC", "logLik" = "logLik",
                            "Pseudo R2" = "pseudo.r.squared"),
             scale = TRUE, robust = TRUE, digits = 2, vifs = TRUE,
             note         = "{stars}. Erros padrões robustos clusterizados",
             title = "Table x - PCAOB treatment",
             to.file = "docx",
             file.name = "PCAOB treatment.docx")
```

# Small Profits (SP)

Criar bases das variáveis dependentes

```{r}
bdSP <- Base_DiD %>%
  dplyr::filter(SP != "NA")
```

## PSM para SP

### Estimando as empresas comparáveis

Identificando as empresas da amostra com maior probabilidade de serem comparáveis

```{r}
PSM_SPPais <- glm(Pais ~ SIZE_w1 + ROA_w1 + LEVER_w1 + BTM_w1 + CFO_control_w1, family = binomial(), 
               data = bdSP)

summary(PSM_SPPais)
```

Estimando as empresas comparáveis

```{r}
Predict_PSM_SP <- data.frame(pr_score = predict(PSM_SPPais, type = "response"),
                     "EUA" = PSM_SPPais$model$Pais)
head(Predict_PSM_SP)
```

Visualizando os grupos previstos

```{r}
#Visualizando os grupos previstos em dois quadros
Grafico_Distr_original <- paste("Distribuição dos grupos:", c("EUA", "JAPAO"))
Predict_PSM_SP %>%
  dplyr::mutate(EUA = ifelse(EUA == 1, 1, 0)) %>%
  ggplot(aes(x = pr_score)) +
  geom_histogram(color = "white") +
  facet_wrap(~EUA) +
  xlab("Score PSM") +
  theme_bw()

#Visualizando os grupos previstos no mesmo quadro
Grafico_Distr_original <- paste("Distribuição dos grupos:", c("EUA", "JAPAO"))
Predict_PSM_SP %>%
  dplyr::mutate(EUA = ifelse(EUA == 1, "EUA", "JAPAO")) %>%
  ggplot(aes(x = pr_score, fill = EUA)) +  # Use 'fill' para cores diferentes
  geom_histogram(color = "white") +
  xlab("Score PSM") +
  theme_bw() +
  labs(title = Grafico_Distr_original) +
  scale_fill_manual(values = c("EUA" = "blue", "JAPAO" = "red"))  # Defina cores personalizadas
```

### Usando a função "matchit" para criar a nova base

```{r}
mod_PSM_SPPais <- matchit(Pais ~ SIZE_w1 + ROA_w1 + LEVER_w1 + BTM_w1 + CFO_control_w1,
  method = "nearest", data = bdSP)

bdSPPais <- match.data(mod_PSM_SPPais)

dim(bdSPPais)
```

## Regressão para SP

Criando base com as variáveis que iremos utilizar e sem NAs

```{r}
bdSPPais <- bdSPPais %>%
  dplyr::select(Ticker, Year, NAICS_Sector, Country, SP, Evento, Pais, SIZE_w1, ROA_w1, LEVER_w1, BTM_w1, CFO_control_w1, Big4, AUD, LCS, COV) %>%
  filter_all(all_vars(!is.infinite(.))) %>%
  na.omit()

summary(bdSPPais)
```

### Estimando o modelo (LOGIT)

```{r}
mod_SPPais_final <- pglm(factor(SP) ~ factor(Pais)*factor(Evento) + SIZE_w1 + LEVER_w1 + BTM_w1 + CFO_control_w1 + factor(AUD) + factor(LCS) + factor(COV),
                   family = binomial('logit'),
                   model = "pooling", 
                   data = bdSPPais, 
                   index = c("Ticker", "Year"),
                   start=NULL, 
                   method = "bfgs", 
                   print.level = 3, 
                   R = 5)

mod_SPPais <- glm(factor(SP) ~ factor(Pais)*factor(Evento) + SIZE_w1 + LEVER_w1 + BTM_w1 + CFO_control_w1 + factor(AUD) + factor(LCS) + factor(COV), 
              family = binomial('logit'), 
              data = bdSPPais)

bdSPPais$probSPPais <- mod_SPPais$fitted.values

summary(mod_SPPais)

summary(mod_SPPais_final)
```

Comparando os coeficientes entre os modelos

```{r}
export_summs(mod_SPPais,mod_SPPais_final)
```


### Estimando a Razão de Chances

Exemplo: para cada variação unitária (aumento) em ETR_log_w1, diminui em 0,81% ((1,941-1)*100) as chances da ocorrência de a empresa ter uma pontuação ESG.
```{r}
logitor(factor(SP) ~ factor(Pais)*factor(Evento) + SIZE_w1 + LEVER_w1 + BTM_w1 + CFO_control_w1 + factor(AUD) + factor(LCS) + factor(COV), 
        data = bdSPPais)
```

### Determinando o Intervalo de Confiança

*A determinação do intervalo de confiança do modelo proposto é relevante para que seja analizada a estimativa do intervalo de predição do coeficiente da variável independente, a um nível de confiança de 95%. Desta forma, em 95% dos casos, o parâmetro dos coeficientes estará dentro deste intervalo.
```{r}
exp(cbind(OR=coef(mod_SPPais_final), confint(mod_SPPais_final)))
```

### Analisando a Curva ROC e AUC(*Area Under Curve*)

```{r}
require(pROC)

roc_SPPais=plot.roc(bdSPPais$SP,fitted(mod_SPPais))
```

```{r}
plot(roc_SPPais,
     print.auc=TRUE, 
     auc.polygon=TRUE, 
     grud=c(0.1,0.2),
     grid.col=c("green","red"), 
     max.auc.polygon=TRUE, 
     auc.polygon.col="lightgreen", 
     print.thres=TRUE)
```

### Obtendo a matriz de confusão

Specificidade e Sensitividade

```{r}
bdSPPais$pdata <- as.factor(ifelse(predict(mod_SPPais, newdata = bdSPPais, type = "response")>0.714,"1","0"))

confusionMatrix(bdSPPais$pdata, factor(bdSPPais$SP), positive="1")
```

### Exportando os resultados com erros padrões robustos clusterizados

```{r}
export_summs(mod_DSSPais_EA, mod_SPPais_final, scala = F,
             model.names = c("EQ.1-DSS","EQ.1-SP"),
             error_pos = c("same", "below", "right"),
             bold_signif = 0.1,
             borders = 2,
             outer_borders = 2,
             statistics = c(N = "nobs.1", R2 = "r.squared",
                            adj.R2 = "adj.r.squared", p.value = "p.value", 
                            "GL" = "df", AIC = "AIC", "logLik" = "logLik",
                            "Pseudo R2" = "pseudo.r.squared"),
             scale = TRUE, robust = TRUE, digits = 2, vifs = TRUE,
             note         = "{stars}. Erros padrões robustos clusterizados",
             title = "Table x - PCAOB treatment",
             to.file = "docx",
             file.name = "PCAOB treatment.docx")
```

